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Abstract: We consider an infinitely large population under stabilising selec¬ 
tion and mutation in which the allelic effects determining a polygenic trait 
vary between loci. We obtain analytical expressions for the stationary ge¬ 
netic variance as a function of the distribution of effects, mutation rate and 
selection coefficient. We also study the dynamics of the allele frequencies, fo¬ 
cussing on short-term evolution of the phenotypic mean as it approaches the 
optimum after an environmental change. We find that when most effects are 
small, the genetic variance does not change appreciably during adaptation, 
and the time until the phenotypic mean reaches the optimum is short if the 
number of loci is large. However, when most effects are large, the change of 
the variance during the adaptive process cannot be neglected. In this case, 
the short-term dynamics may be described by those of a few loci of large 
effect. Our results may be used to understand polygenic selection driving 
rapid adaptation. 
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In the study of fast adaptation, polygenic selection may be more im¬ 
portant than selection on single genes. At single genes, strong selection 
driving fast adaptation generally leads to the rapid fixation of beneficial 
alleles or at least to huge allele frequency shifts between populations. Clas¬ 


sical examples of this type o 


i n moths (van’t Hof et al. 


(Daborn et al. 


adap tation are the case of industrial melanism 


2 0111 1 and insecticide resistance in Drosophila 


20021). However, in many (if not most) cases of fast adap¬ 


tation, such as in island lizards that are able to adapt very quickly to a 
changing; veg etation (e.g. revealed in an experimental evolution study by 


Kolbe et al 


(120121) ). small shifts of allele frequencies at many loci may be 


sufficient to move a phenotype toward a new optimum under changed envi¬ 
ronmental conditions. 

There is a large and growing body of literature on the detection of adap- 
tive signatures in molecul ar po pulation genetics. Following pioneering work 


of 


Maynard Smith and HaighI (1974), the impact of positive selection on 


neutral DNA variability (selective sweeps) has attracted much interest. This 
theory has been applied to huge datasets that emerge from modern high- 
throughput sequencing. A large number of statistical tests have been de¬ 
veloped t o detect sweep signa l s and estimate the frequency and strength o 


selec 


ion (K im and Steph an 


2002 


Nielsen et al. 


2005 


Pavlidis et al. 


2010). However, most theory so far excludes the phenotypic side of the 


adaptive process (except for fitness). Usually, selection is simply modeled as 
a constant force that acts on a new allele at a single locus. This is in strik- 
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ing contrast to the classical phenotype-ba sed models of adaptat ion that are 


successfully used in quantitative genetics (BARTON and KE IGHTLE Y. 


2002 ). 


These models typically assume that adaptations are based on allele frequency 
shifts of small or moderate size at a large number of loci. Also, adaptation 
does not require new mutations, at least in the short term. Instead, selection 
uses alleles that are found in the standing genetic variation. Genome-wide 
data of the past few years show that this quantitative genetic view is relevant. 
In particular, association studies confirm that quantitative traits are typically 
highly polygenic. High heritabilities most probably result from standing ge¬ 
netic variation at a large number of loci with small individual effect. Also, 


local adaptation to environmenta 


at multiple loci (Hancock et al. 


dines involves moderate frequency shifts 


201 0). As a consequence, there is grow¬ 


ing evidence that the molecular scenario of sweeps only covers part of the 
adaptive process and needs to be revised to include polygenic selection. 

As genome-wide association studies (GWAS) yield information about the 
distribution of single nucleotide polym orphisms (SNPs) relevant to quan¬ 


titative traits ( VlSSCHER. et al 


201 21 ). it is important to understand the 


models of polygenic selection in terms of the frequency changes of molec¬ 
ular variants, i.e. in terms of population genetics. So far, however, the 
dynamics of only very simple poly genic models have been studied and ap¬ 


plied to data (e.g. 


Turchin et al. 


the dynamics of a quantitative 
that was originally proposed by 


(2012)). In this article, we analyze 
rai t giv en by_a much more general model 


Wrig htI (119351) and recently re-visited by 
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de Vladar and Barton 


(2014). These authors consider an infinitely large 


population evolving under stabilizing selection and mutation. Following the 
observa ti ons from many empirical studies (particularly from biomedicine, see 


VlSSCHER et al. 


(2012)), they assume that the effects are locus-dependent. 


Furthermore, they consider additivity of the effects and linkage equilibrium 
between loci. 

Here we study this model and obtain some analytical results on the sta¬ 
tionary genetic variance and also the dynamics of the phenotypic mean. We 
show that the stationary genetic variance may exhibit nonmonotonic depen¬ 
dence on the shape of the distribution of effects. We also study how the 
trajectories of the allele frequencies and the mean trait respond to a sudden 
environmental shift. When most effects are small, as is th e case in expe rt 


ments on Drosophila (MACKAY^ 


2001 


Goddard and Hayes 


2004 1. in livestock ( Hayes and Goddard 


2009) and for human height (VlSSCHER . 


2008j). 


a simple analysis shows that the magnitude of the deviation of the pheno¬ 
typic mean from the optimum decays roughly exponentially with time and 
approaches zero over a time scale that is inversely proportional to the initial 
genetic variance. When most effects are large, the short-term dynamics of 
the mean and variance can be understood by considering a few loci with large 
effects. 


MODEL WITH LOCUS-DEPENDENT EFFECTS 


We consider the t-locus model recently analysed by 


de Vladar and Barton 
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(1201 4) where each locus is biallclic. The + allele at site i has frequency p t 
while the — allele occurs with frequency q t — 1 — pi. The effects are assumed 
to be additive so that the trait value is z = Yh\=i sgn(i) lu where sgn(i) = ±1 
denotes the sign of the genotypic value of locus i and 7 * > 0 is the effect of 
the allele at the it\i locus. The loci are assumed to be in linkage equilibrium 
so that the joint distribution of effects at the loci factorises. As a result, 
the nth cumulant c n of the phenotypic effect, obtained on averaging over the 
population distribution, can be written as the sum over the corresponding 
quantities at individual loci. The first t hree cumu lants viz. mean ci, variance 


C 2 and skewness C 3 are given by (IBurgerI . 


199l|) 


Cl = 


c 2 = 


c 3 = 


~ Qi) 

i=l 

£ 

2 ^ mi 


(ia) 

(lb) 

(lc) 


1=1 


The allele 


by (Barton, 


frequ ency evolves in time under selective pressure and is given 


19861 ) 


dpi 

dt 


~Pi(t + 1) 


Pi(t) 


PiQi dw 
2w dpt 


( 2 ) 


where w is the average fitness of the population. For large £, as the trait 


6 















value z of an individual can be treated as a continuous variable, from (Hal) 
and fllbl) . we obtain 


w — I dz p(z) w(z) = 1 -(c 2 + (Aci) 2 ) « e 2 

./—^ 2 


-f ( C2 + (A C1 ) 2 ) 


(3) 


where the approximate equality sign holds because s is assumed to be small. 
In the above expression, w(z) = 1 — (s/2)(z — z a ) 2 is the fitness distribution of 
the phenotypic trait under stabilising selection, z Q the phenotypic optimum 
and Aci = c\ — z a the mean deviation from z 0 - Thus the maximum fitness 
(namely, one) is obtained when the population is at the phenotypic optimum 
and has no genetic variance. Inserting equations © and (| 3 ]) in (j2j) and 
accounting for mutations, we obtain the following basic equation for the 
evolution of allele frequencies, 

'-jjj- = PiQi + Pi) + ~Pi) , i = , (4) 


where /i is the probability of (symmetric) mutation between t he + and — 


allele at locus i. Note that the equation (1) of 


deVladar and Barton 


12 0141 1 is obtained by replacing s by 2 s in the above equation. 


On the right hand side (RHS) of (141) , the first term (in the first parenthe¬ 
sis) expressing the mean deviation from the optimum corresponds to direc¬ 
tional selection toward the phenotypic optimum: if the mean is above (below) 
the optimum, the allele frequencies decrease (increase). However, once the 
phenotypic mean is sufficiently close to the optimum, stabilising selection 
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(described by the second term) takes over. 

One of the difficulties in solving (j3J) is that it involves the mean c\ which 
depends on all the allele frequencies. Moreover, it has been shown that the 
differential equations for the cumulants do not close: each one not only in¬ 
volves two higher cumulants, but also contains terms that can not be writte n 


in terms of other cumulants (IBarton and TURELLI 


19871 


B urgerI . 


199 lh . 


GENETIC VARIANCE IN THE STATIONARY STATE 


In the stationary state in which the left hand side of flj]) vanishes, if the 
mean c* = z 0 , the allele frequency p* has three solutions, namely 1/2 and 
(1 ± y/l “ (7/7;) 2 ) A where 7 = 2 \/2p/s. The latter two solutions are 
stable for 7 \ > 7 , and therefore the allele frequency is close to fixation when 
the effects are large. For 7 \ < 7 , the effects are small and the stationary 
state so lution p*_ _ 1/2 is t he on ly stable solution for the allele frequency 


f DE Vl adar and Ba rtonI . 


2 0141 1. From th ese resu lt s, t h e sta tionary genetic 


variance (llbh is easily seen to be (IDE Vladar and Barton 


2014]) 


* 4 /^ 1 

C * = T m+ 2 


Et?. 


(5) 


7i<7 


where, for large £, the number of effects larger than 7 is given by ni = 
£ /.°° d'f p(y) with p(y) being the distribution of effects. Thus the genetic 

















variance in the stationary state can be neatly written as 


c 9 = - 
2 2 


7 2 / d/y p(yy) + / d'y y 2 p(y) 


( 6 ) 


where the first (second) term is the contribution from loci with large (small) 
effects. 

The gamma distribution p(yy) ~ 7 fe_ 1 e _fe7 /^ with shape parameter k > 
0 and mean 7 has been emp loyed to fit the distribution of QTL effects 


(Hayes and Goddard, 


2001 ). This distribution is L-shaped for k < 1 and 


bell-shaped for k > 1, while for k — 1, it is an exponential function. For the 
gamma distribution, the stationary genetic variance (J 6 ]) is given by 


Co = 


f?7 2 


T(fc, k^f r ) T (2 + k) — T (2 + k, k'jr) 


m 


7 2 r k 2 T(k) 


(7) 


where 7 r = 7/7 and T(a,b) = dt i a_ 1 e _t is the incomplete gamma func¬ 
tion. For the special case of exponentially distributed effects {k = 1 ), (J7J) 
simplifies to give c* 2 = ^ 7 2 (1 — e _7r (l + y r )). 


For y > 7, 


4 pi/s (Tu relli 


le House of Cards (HoC) variance, namely c 2 = ^y 2 /2 = 


19841 ). is obtained when 7 is finite but 7 —> 00 for all fc. In 


the opposite case (7 7 ), we have c 2 = £ 7 2 (fc + l)/(2/c) which depends on 

the shape of the distribution of effects. For fixed mean 7 , the genetic variance 
increases monotonically with the scale 7 since a larger mutation probability 
increases the variance. If instead the distribution mean 7 is increased keeping 
7 fixed, the variance increases with 7 toward the HoC value. This is because 
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for fixed k, the width of the distribution increases with 7 and therefore larger 
effects can be accessed. 


The stationary genetic variance has been computed nume rically f or vari ¬ 


ous s hape parameters when 7 = 0.063 ,7 = 0.1 and £ = 1000 in 


de Vladar and Barton 


(120141 1. From QT|), the variance for k = 1,2,10 and 100 is found to be 
1.32,1.57,1.94 and 2, respectively, which agrees well with the numerical data 
in their Figure 5. To understand how the variance depends on the shape pa¬ 
rameter, we first note that with increasing k (and fixed 7 ), the width of the 
gamma distribution decreases. For large k, if 7 > 7 , the variance saturates 
to the HoC variance since almost all loci have large effects with narrow dis¬ 
tributions while in the opposite case, most effects are small and the variance 
tends to £f 2 (see Fig. [[]). For small k, irrespective of whether 7 is above 
or below 7 , we find that most effects are small. To see this, consider the 
fraction f s = 1 — ( ni /£) of loci with small effects which is given by 


_ (HO* 

Is 


(k-l)\Jo 


dx x k -'e- kr/rX = 1 - T ( k ’ k ^ r 


m 


( 8 ) 


If k < ( 7 r ) _1 , the above equation yields f s ~ {k^ r ) k /k\. Then for finite y r 
when k —> 0, we find that n s —>• £ for any 7 r , as claimed above. To summarise, 
as shown in Fig. Q] for 7 < 7 , the variance increases with k toward the HoC 
variance, while for 7 > 7 , both c 2 and n s are nonmonotonic functions of k. 

When the effects are chosen from an exponential distribution, the fraction 
f s — 1 — e _7r . On eliminating 7 r in favor of f s in 03) for k = 1, we find the 
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relative contribution of loci with small effects to the total variance to be 



2/s + (1 — f s ) ln(l — / s )(2 — ln(l — f s )) 
2fs + 2(1 — f s ) ln(l — f s ) 


( 9 ) 


which increases as f s /3 for small f s and approaches unity as f s increases 
toward one. The above expression shows that if 10% of the effects are small, 
their contribution to the variance is merely 3%, which increases to 21% when 
f s is one half. To obtain an equal contribution from small and large effects, 
a disproportionately large fraction (~ 83%) of small effects is required. We 
are unable to obtain an analytical expression analogous to (jU]) for arbitrary k 
since f s is not a simple function of y r (see (J8J) above). However, a numerical 
analysis using (J7J) and (JHJ) shows that for the same value of f s , small effects 
contribute more to the total genetic variance as the distribution of effects 
gets narrower. For f s = 0.1, the relative contribution is found to be 2%, 3% 
and 5% for k = 1/2,1 and 2, respectively. To obtain an equal contribution 
from loci with small and large effects, f s = 0.89,0.83 and 0.77 is needed for 
the shape parameters k = 1/2,1 and 2, respectively. 


DYNAMICS OF THE ALLELE FREQUENCY 


We now turn to a description of the allele frequency dynamics, and will 
consider the situation when the phenotypic optimum is suddenly shifted. As 
mentioned earlier, due to the term Aci on the RHS of (|4lh all the frequencies 
are coupled which makes it hard to obtain an exact analytical solution of the 
allele frequency dynamics. However, under certain conditions, it is a good 
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approximation to consider only the C\ term in (j3J) for the initial dynamics 
and the rest of the terms for long-term evolution. 

To see this, we first note that since 0 < Pi < 1, the mean |ci(t)| < JA 7 * ~ 
7 't. For independent and uniformly distributed initial frequencies, as the 
average initial frequency is one half, the leading order contribution (in £) to 
the initial mean is zero. The initial variance is, however, nonzero which gives 
the typical initial mean |ci(0)| ~ 7 \/£. When the phenotypic optimum z Q < 
and the number of loci is large, the initial value |Aci( 0 )|/ 7 i > 1 . 

Thus at short times, we can neglect \2pi~ 1| (which is bounded above by one) 
and the mutation term in comparison to the term in (|4]). At large 

enough crossover time t x , as explained below, the mean deviation is close to 
zero and the reverse condition holds; i.e. 2 |Aci(£)|/ 7 j ^ |2 Pi — 1| in (j3J), and 
we may set Aci ~ 0 for later evolution. Biologically these considerations 
mean that initially the effects are weaker than the mean trait deviation, but 
as the population adapts due to directional selection, the deviation of the 
mean from the phenotypic optimum becomes smaller than the effects. 

The above argument applies not only to uniformly distributed initial fre¬ 
quencies but in more general settings as well where |ci( 0 )| ~ 7 1 by replacing 
\fl by i. Here we will focus on the dynamics of the allele frequency when 
the optimum is suddenly shifted to a new value Zf(< £ 7 ), starting from the 
population which is equilibrated to a phenotypic optimum value z 0 . In this 
situation, as the initial frequency is close to one half when 7 \ < 7 , the fre¬ 
quency \2pi(0) — 1 | is obviously negligible compared to Aci( 0 )/ 7 j, whereas 
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for 7 i > 7 , | 2 pi( 0 ) 


]J i s c l ose to one s ince the initial frequency is close to 


either zero or one ( PE VLA DAR and Ba rton 


20141). 


When most effects are small: The effects at most of the loci can be 
smaller than the scale 7 either if k is large and 7 < 7 , or if k is small. 
Then for most loci, at short times, the full model defined by (J3J) can be 
approximated by 


dpi A 

— = -s^iPiqi/Xci 

at 

dCr, 


dt 


= -sAcic n+ i , n > 1 , 


(10a) 

( 10 b) 


where the last equation for cumulants is obtained from the results of 


Burger 


( 199ll ). Equation (llObl) for n = 1 shows that the magnitude of the mean de¬ 
viation decreases with time, and for the phenotypic optimum smaller than 
the maximum attainable value of the mean ( z 0 <C 7 T), the trait mean be¬ 
comes close to the phenotypic optimum at large times (see Fig. [ 2 ja) . We now 
assume that the variance c? is independen t of time and stays at its initial 


value c 2 (0) ( CHEVIN and HOSPITAL . 


2008). As explained in Appendix lAl 


this approximation is good when a combination of the initial cumulants is 
small (see also Fig. [2b). This allows us to solve (HOaj) and (llObj) . and we 
immediately find that 


Aci(t) 

Pi{t) 


Aci( 0 )e- C2(0)st 

Pi(0) 


, s , s 7 ,;Aci(0) , 1 

Pi(0) + qi(0)e C 2 (Q) 


( 11 ) 

( 12 ) 
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Equation (fTTD shows that the mean deviation approaches zero over a time 
scale t x (sc 2 ( 0 )) \ 

Next we analyse the long-term evolution of the allele frequencies. As 
Fig. [21a, shows, there is a small but nonzero mean deviation AcJ in the sta¬ 
tionary state. Taking this into consideration and accounting for the other 
terms in TP) , for t > t Y . we can write 

% = - 2 Pi + —") + M 1 ~ 2 Pi) ■ ( 13 ) 

ot 2 7 i 

For AcJ = 0, the above equation can be easily solved to give 

(14) 

where 


Mi(t) 


4 Pi(t x ) - 4Pi(t x ) + rrii f 2 
( 2 Pi(t x ) - l ) 2 


5 t ^ 1 


(15) 


m i = ( 7/7 i) 2 and Pi(t x ) is obtained from (TT2l) . We check that the stationary 
state solutions ( 1±\/1 — rrii )/2 and 1/2 are obtained from the above result for 
TYii < 1 and > 1, respectively. Furthermore, the solution p[ + \t) is obtained 
for pi(t x ) > 1/2 and p['\t) for p l {t y ) < 1 / 2 . 

Figures [2] and [3] show a comparison between the numerical solution of 
(J4j) and the approximation described above, when the initial condition is the 
stationary state of the population equilibrated to a phenotypic optimum z Q . 
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The initial mean deviation Aci(O) is seen to be close to — Zf, and the initial 
variance 02 ( 0 ) = c\ for the zero mean deviation is 0.0967, which is close 
to the value 0.131 obtained from the set of effects used in Figs. E] and El 
As Fig. [2] shows, the dynamics of the mean deviation are captured well by 
(TITD and approach a stationary value close to zero (Ac^ ~ —0.016) in about 
1500 generations. The variance also evolves with time, but the change is not 
substantial and the approximation c^it) ~ C 2 ( 0 ) is good over the time scale 
directional selection toward the phenotypic optimum operates. Equation (TITjl 
also indicates that directional selection toward the optimum will occur faster 
when the initial variance is large since t x ~ 1 jc* 2 . 

Since the stationary genetic variance displays a nonmonotonic depen¬ 
dence on the shape parameter k of the gamma distribution (see Fig. [[]), the 
relaxation time for the mean deviation is expected to decrease and then in¬ 
crease with increasing k. Indeed, as the inset of Fig. Ek shows, the difference 
ci(i) — cl (which, by definition, is zero in the equilibrium state for all k) 
equals a reference value —0.05 at time 360, 340 and 370 for k — 1,5 and 20, 
respectively. 

The allele frequency dynamics are shown in Fig. [3j We see that while 
the short-term dynamics can be accurately described by (fT2l) for loci with 
effects smaller than or close to the distribution mean, there is a substantial 
difference when the effects are larger than the mean. This is because for such 
loci, the initial frequency is not close to half and the term involving g* — pi 
on the RHS of (j3J) can not be neglected. For t > t x , the long-term behavior 
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described by is shown with Ac* = 0 and the actual mean deviation. 
When most effects are large: When 7 > 7 and k is large, the number 
of loci with large effects is also large, and the initial allele frequencies are 
close to either zero or one. In this parameter regime, both the variance and 
the skewness may change appreciably during directional selection toward the 
optimum, and the constant-variance approximation discussed above is not 
suitable. However, at very short times when Aci(f) is close to its initial 
value, the solution (fl 2 | ) for the allele frequency gives 


Pi(t) ~ 


1 


1 + 


9i(0) Aei (0W 

vm e 


(16) 


From the above equation, we first note that the allele frequency at large- 
effect loci changes fast as expected intuitively. Equation (flfili also shows that 
for Aci(O) < 0, the allele frequency quickly increases toward unity, if the 
initial frequency is close to unity and therefore does not contribute to the 
dynamics of the variance or skewness. Thus to understand the short-term 
dynamics, we need to focus our attention on large-effect loci with low initial 
allele frequency for negative initial mean deviation. Similar remarks apply 
to the situation when Acq (0) is positive where the large-effect loci with high 
initial frequency determine the dynamics. 

In the following, we assume that Aci(0) < 0 and consider the time evo¬ 
lution of the allele frequency P of the largest effect locus with lowest initial 
frequency. Figure [I] shows that the allele frequency P sweeps to fixation, but 
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the frequency of the next relevant locus (i.e. the next largest effect locus 
with low initial frequency) does not. In such a case, we can approximate the 
mean c± and variance C 2 by the contribution from the frequency P with effect 
r, and obtain 


Cl(t) 

» 2T(P-P 0 ) + Ci(0) 

(17a) 

C 2 [t) 

~ 2T 2 (PQ - P 0 Qo) + c 2 (0) , 

(17b) 


where P 0 = P( 0). Then using the above expression for the phenotypic mean 
in (jlD and neglecting mutations (since most effects are large), we get 


dP 

~dt 


—sT 2 P(l — P)(P + a) , 


(18) 


where 


a = 


T + 2Aci(0) 

2f 


-2P 0 • 


(19) 


We thus find that the allele frequency P is a solution of the following equation 


(P/P 0 ) 1+a _ sr 2 a(l+a)t P + a /'Of)') 

(Q/Qo) Q Po + a • 1 J 

An explicit solution of (12(1 seems hard to obtain since a is in general not an 

integer. However, for large and negative a, the above equation yields 


P 


1 


1 + 


Qo_ p — sr 2 |a|t 
Po ^ 


( 21 ) 
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Thus for Zf T, the frequency P sweeps to fixation in a time of order 
(sTz f )~\ 

Figure 3] shows the allele frequency of the largest effect locus with low¬ 
est initial frequency obtained using (J3J). It agrees reasonably well with the 
solution of (l20p and the expression (l2Tj) where a = —1.43. In Fig. [51 the 
dynamics of the first two cumulants given by (Ilap and (llbH are compared 
with the approximate expressions (I17ap and (117bf) . respectively, and we see 
a good agreement. 

A detailed numerical analysis of the set of parameter values of Fig. [4] 
suggests that the dynamics of this example can be understood by considering 
one, two or three of the largest-effect loci. When only the largest-effect locus 
was required, the second largest effect was much lower than the largest one. 

DISCUSSION 


One of the fundamental questions in adaptation is whether the adap¬ 


tive pr ocess 


effect ( OrrI. 


is governed by many loci of small effect or few loci of large 


200511 . However, which effects are small, and which large? 


deVladar and Barton 


(120141 1 have provided a scale 7 ~ \fjjjs for the 


size of effects, which is a function of basic population genetic parameters, 
namely mutation probability p and selection coefficient s, relative to which 
an effect is defined as large or small. For a given distribution of effects, as¬ 
sumed here to be a gamma distribution with mean 7 and shape parameter 
k, an effect is small (large) if it is below (above) 7 . But for fixed 7 and 7 , 
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whether most or a few effects are small depends on the shape parameter k: 
for small k, most effects are small, but for large k, the number of small effects 
depends on the ratio 7 r = 7 / 7 . 

Genetic variance in stationary state: Here we have provided analytical 
expressions for the stationary genetic variance c\ when the effects are locus- 
dependent. We find that when most effects are small, is a nonmonotonic 
function of the shape parameter of the gamma distribution (going through a 
maximum for intermediate values of k ; see Fig. [[]). In contrast, it increases 
monotonically when most effects are large. As Fig. [T] shows, when the shape 
distribution is narrow, large (small) effects contribute most to the variance 
when 7 r > 1(< 1). However for broad distributions, although the number of 
small effects is large, small effects do not contribute much when y r < 1. The 
HoC variance is obtained irrespective of k when y r —> 0 si nce all loci have 


l arge effect in this limit. As noted previously (IDE Vlad ar and Ba rton 


2 0141 ). HoC provides an upper bound on the genetic variance. 


Dynamics when most effects are small: As the distri 


Dution of QTLs_ 


measured in experimental anc 

natura 

populations 

(Mac kay 

2004; 

Hayes and Goddard 

2001 ; 

Goddard and Hayes, 

2009; 

VlSSCHER 

, 

2008 

) find most effects to 


be small, it is important to study this situation in detail. Here we have 
obtained analytical expressions for the dynamics by assuming the genetic 
variance to be constant. Although the fact that the variance does not 
change much in ti me when most e ffects are small was observed numerically in 


de Vladar and Barton 


(20141). an explanation of this behavior was not 
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provided. Here, as explained in Appendix 1X1 it is a good approximation 
to assume the variance to be time-independent provided the product of the 
initial values of the mean deviation and skewness is small. 


In the absence of mutations, 


C hevin a nd Hospital! (2008|) have con¬ 


sidered the effect of background with a time-independent genetic variance 
on the frequency at a single focal locus. Their results match the ones 
obtained here using the short-term dynamics model with directional selec¬ 
tion _only; in partic ular. ([Till and (1T2]) match the results ( 21 ) and (25) of 


Chevin a nd H ospital! (120081 ). respectively, on identifying their parameters 


cu 2 and a with 1 /s and 7 from this study. 

Our basic result concerning the dynamics of the phenotypic mean is that it 
relaxes over a time scale which is inversely proportional to the initial variance. 
Since the variance is of order £, we thus have the important result that the 
mean approaches the optimum faster if a larger number of loci is involved. 
Moreover, this time depends nonmonotonically on the shape parameter of 
the gamma distribution. Note that the phenotypic mean deviation relaxes 
to zero when the phenotypic optimum is far below the upper bound 7 1 on 
the phenotypic mean. However, when the phenotypic optimum exceeds the 
maximum typical value of the mean, such that the mean deviation remains 
substantially different from zero at late times, (llObl) shows that all higher 
cumulants vanish at the end of the phase of directional selection. 

Dynamics when most effects are large: When the initial mean deviation 
is moderately large (and negative), the genetic variance changes by a large 
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amount over the time scale directional selection occurs and the dynamics 
can be understood by considering a few loci whose effect is large but initial 
frequency is low. 

However, for larger mean deviations (but smaller than £ 7 ), a few large- 
effect loci do not completely capture the dynamics of the mean and the 
variance. As Fig. ISII shows, the initial increase of the absolute mean deviation 
and the transient rise of the variance can be explained by considering the 
large-effect locus. At later times, however, as the change in variance is small, 
we can use the constant-variance approximation to understand the dynamics 
of the phenotypic mean deviation until it nearly vanishes. The constant- 
variance approximation can also be used when the initial mean deviation is 
sufficiently small (see Fig. IS 2 |) . 

Applications: The approximations presented here hold for the short-term 
evolution of phenotypic traits and allele frequencies. This means that our 
results may be used to understand polygenic selection driving rapid adap¬ 
tation. In this respect, our most important result is that the mean of a 
phenotypic trait may respond faster to a sudden environmental change when 
the number of loci is large and most effects are small. 

Evidence for rapid phenotypic evolution has been reported in recent years 
from several groups of organisms. For instance in Drosophila subobscura , 


latitudinal dines of wing size h ave been 


species colonized America (Huey et al. 


brined within 20 years since this 


200QI). Similarly, in held experi¬ 


ments in which lizard populations were newly established on small islands in 
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the Bahamas, t 


le hindlimbs a dapte d very quickly to the different vegetations 


on the islands (Ko lbe et al. 


2 012 ). To our knowledge, however, data from 


GWAS are not yet available in these cases. 

The theory presented here can also be applied to the large amounts of 
GWAS data that have been gathered in model species such as humans and 
Drosophila. To analyse the observed allele frequency shifts in_SNPs associ¬ 


ated with quantitative traits, such as h uman heig 


and cold tolerance in Drosophila (IHuang et al. 


it 


Turchin et al. 


20121 ) 


in this study provide a more general th eoret ical basis t 
equations used in previous analyses (e.g. T urchin et al . 


20121 ). the results derived 
lan the dynamical 


( 2012 )) 


Open questions: The analytical calculations in this article work when the 
phenotypic mean at the equilibrium coincides exactly with the optimum. 
However, in the stationary state, there is a small but nonzero mean de¬ 
viation due to which the long-term dynamics are not accurately captured, 
especially when effects are small, as shown in Fig. [3)3. An improved calcu¬ 
lation of the dynamics that takes a nonzero mean deviation into account is 


ations (leading to sel< 

active sweeps) in this model ( 

Chevin and Hospital. 

2008; 

Pavltdts et al. 

2012; 

WOLLSTEIN and Stephan, 

2014 

)• 


Another open question concerns the generality of our results presented 
here. The current model perhaps oversimplifies biological reality in that it 
neglects genetic drift and assumes additive effects, symmetric mutations and 
free recombination between loci. It can be shown that the current model 
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(neglecting mutation) can be derived from the classical symmetric viabil¬ 
ity model with arbitrary position of the optimum under the quasi-linkage 
equilibrium assumption. In the latter model the probability of selective hxa- 


tion h as bee n studied numerically 


Wollstein and Stepha n. 


for up to eight loci (IPavlidis et al. 


2012 


20141 ). However, at present we are lacking an 


analytical understanding of the role of recombination in this model and how 
it relates to the high-recombination limit represented by our current model. 
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APPENDIX 


A Validity of constant-variance approxima¬ 


tion 


As explained in the main text, we obtain (TTTT) when the variance is assumed 
to be constant in time. To see under what conditions this approximation 
holds, we find a correction to the variance by assuming that the skewness C 3 
is nonzero and time-independent (see the inset of Fig. Od). Using (ITTj) on the 
RHS of (jlObj) for n = 2, we immediately get 


c 2 (t) = c 2 (0) [l-C(l-e^ st )} 


(A.l) 


where C = Aci(0)c 3 (0)/c|(0). Plugging the above solution into f)10b[) for 
n — 1 gives the mean deviation as 


Aci(i) = Aci(O) exp [—sic 2 (0)(l — C) — C (1 — e stC2 ^)] 


(A. 2 ) 


The solution (ITTj) is recovered if the constant C, which depends on the initial 
value of the first three cumulants, is negligible. 

If we start with the initial condition in which the population is equili¬ 
brated to an optimum and most effects are small, since most initial allele 
frequencies are close to one half, the initial variance is substantial and the 
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skewness is close to zero. In this case, the constant-variance approximation 
is expected to work well. However, if most effects are large, since most initial 
allele frequencies are close to fixation, although the skewness remains small, 
the variance also becomes small, thus leading to an increase in C. Then for 
sufficiently small mean deviations, we may also employ the const ant-variance 
approximation when most effects are large. 
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Variance 



ShapeParameter 


Figure 1: Genetic variance in the stationary state as a function of the shape 
parameter k when the effects are distributed according to the gamma func¬ 
tion. The plot shows the total genetic variance (solid), variance due to small 
effects (small dashes), large effects (large dashes), and the fraction of small 
effects (dotted) for (a) 7 = 0.04,7 26 0.08 and (b) 7 = 0.1,7 = 0.05 for 
£ = 1000. The asymptotic values £7 2 (k + l)/(2 k) when 7 > 7 and £^ 2 /2 
when 7 < 7 are also shown (top solid curves). 














MeanDeviation 



(a) 


Variance 



Figure 2: Response to change in optimum when most effects are small. The 
plot shows the results for (a) mean deviation Aci(i) and (b) variance 02 (f) 
and skewness 03 (f) obtained using the exact numerical solution of the full 
model (solid) and the short-term dynamics model (large dashes). The dot¬ 
ted curves show the time-dependent solution (1111) for mean and (1A.1D for 
variance. The parameters are l = 50, s = 0.02,/i = 5 x 10~ 5 ,7 = 0.14 > 
7 = 0.05, z 0 = —0.0012, = 0.5, rti = 5. The effects are chosen from an 

exponential distribution, and the parameter C = —0.01 (see Appendix [Al l. 
The inset in the top figure shows the difference Aci (f) — Ac^ as a function of 
time for the full model when the effects are gamma-distributed with shape 
parameter k = 1 (large dashes), 5 (sig^ll dashes) and 20 (dotted). 
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time 


Figure 3: Response to change in optim um when most effects are small. The 
plot shows the allele frequencies for two representative loci with (a) 7 * = 0.252 
and (b) 7 \ = 0.028 for the full model (solid) and short-term dynamics model 
(large dashes). The dotted curves show the time-dependent solution f|T 2 |) for 
t < t x and (TT4|) for t > t x where t x = 1500. The dashed curve for t > t x is 
the solution of (fl~3l) with Ac^ = —0.016. The other parameter values are the 
same as in Fig. [ 2 j 

















Figure 4: Response to change in optimum when most effects are large. The 
plot shows the exact numerical solution of the full model (solid) and the 
equations (l 20 |) (large dashes) and fl 2 Tf (small dashes) for the dynamics of 
the allele frequency P with the largest effect and lowest initial frequency 
(T = 0.776, P 0 = 3.3 x 10 -4 ). The solid curve at the bottom shows the 
numerical solution of the full model for the frequency of the next relevant 
locus with effect size 0.319 and initial frequency 1.9 x 10~ 3 . The parameters 
are i = 20, s = 0.1,/x = 10~ 5 ,7 ~ 0.028 <C 7 = 0.2, z Q = 7.8 x 10 ~ 5 ,Zf = 
1.5, m = 19. 
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MeanDeviation 



Variance 



Figure 5: Response to change in optimum when most effects are large. Solid 
lines show the mean deviation (TTaT) and variance (jlbl) . while the large dashed 
curves show the contribution to these cumulants from the locus with the 
largest effect and lowest initial frequency (r = 0.776, Pq = 3.3 x 10 -4 ). In 
both cases, the exact numerical solution of the full model is used. The 
numerical solution of (TT8|) (small dashes) is also shown. The other parameter 
values are the same as in Fig. [4j 
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MeanDeviation 




Figure SI: Response to change in optimum when most effects are large. 
Solid lines show the mean deviation (flaj) and variance (llbj) . while the large 
dashed curves show the contribution to these cumulants from the locus with 
the largest effect and lowest initial frequency (r = 0.776, Pq = 3.3 x 10~ 4 ). 
In both cases, the exact numerical solution of the full model is used. The 
numerical solution of (TT8li (small dashes) is also shown. The dotted curves 
show (fTTjl and (lA.ljl for t > 100. The final optimum value Zf = 2.5 and the 
other parameter values are the same as in Fig. [4l 
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Figure S2: Response to change in optimum when most effects are large. Solid 
lines show the mean deviation (flail and variance (llbl) . while the two dashed 
curves show the contribution to these cumulants from the first two relevant 
loci with effect 0.77 (large dashes) and 0.34 (small dashes). In both cases, 
the exact numerical solution of the full model is used. The dotted curves 
show (fTTT) and (IA.If) with C = 0.006. The final optimum value Zf = 0.3 and 
the other parameter values are the same as in Fig. [If 


3 

















